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Abstract 

This paper introduces the basic concepts for physics-compatible discretization techniques. The paper gives a clear 
distinction between vectors and forms. Based on the difference between forms and pseudo-forms and the *-operator 
which switches between the two, a dual grid description and a single grid description are presented. The dual grid 
method resembles a staggered finite volume method, whereas the single grid approach shows a strong resemblance 
with a finite element method. Both approaches are compared for the Poisson equation for volume forms. 
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1. INTRODUCTION 



Mimetic methods aim to preserve essential physical/mathematical structures in a discrete setting. Many of such 

S structures are topological, i.e. independent of metric, and involve integral relations. Since integration will play an 

important role and integration of differential forms is a metric-free operation, we will work with differential forms. 
Formally, differentials forms are linear functionals on multi-vectors, but Flanders, ifTTl p.l], refers to them as 'things 
—— , which occur under integral signs'. Such would not be the case if we were to use vectors, because integration of 

vector quantities is a metric operation. The same holds for vector operations; the grad, curl and div are metric- 
dependent operators, whereas the exterior derivative, which plays a similar role for differential forms, is metric-free. 
2T The important difference between vectors and forms will be explicitly addressed in this paper. 

When integrals over ^-dimensional geometric objects are considered, the orientation of these A:-dimensional ob- 
jects need to be taken into account. If we change the orientation of a point, curve, surface or volume, some integral 
values change sign, whereas others do not. For instance the work Wab of a conservative force along a curve y con- 
necting the points A and B is equal to —Wba, i-e. the work of the same force in the opposite direction along the curve. 
So the physical quantity work changes sign when we change the orientation of the curve. Mass, on the other hand, 
which is the integral of mass density over a volume, does not change sign when we change the orientation. 

Therefore, we need to consider two distinct types of differential forms: Those that do not change sign when 
orientation is reversed, the true forms and those that do change sign, the pseudo-forms. The operator which switches 
between forms and pseudo-forms is called the Hodge-* operator. This operator depends explicitly on the metric. 

Integrals and integral relations can be represented without error in terms of duality pairing between chains and 
cochains. The distinction between integrals of true forms and pseudo-forms requires in principle two grids: One on 
which we represent the integral of a true form and the other grid on which we represent the integral of a pseudo-form. 
The formulation obtained by employing two dual grids resembles staggered finite volume methods. 

An alternative way to implement the action of the Hodge-* operator is to make use of an inner product. In this 
approach only one grid is required. The formulation based on a single grid approach leads to a finite element method. 
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In this paper we try to explain this structure in more detail and show in a specific example that the single grid 
approach and the dual grid approach lead to equivalent solutions. By making the clear distinction between topological 
concepts (integrals and discrete integral representations) and metric dependent operations (Hodge-*), it is very easy 
to switch between orthogonal grids and curvilinear grids. It will be shown that single grid and dual grid methods give 
equivalent solutions in a specific problem. 

Commuting relations between the discretization (mimetic projection) and operations at the continuous and discrete 
level will play an important role in order to ensure that the 'discrete system behaves just like the continuous system'. 

Throughout this paper the basic idea will be highlighted by putting statements in a box and the main idea of this 
paper is: 



A discrete representation of a physical system will display the same structure/dynamics when when discrete 
operators and the continuous operators commute with the projection of the infinite dimensional space onto 
the discrete space. 



The same idea has been put forward in many different papers. Early work in this field was reported by Branin, |[8). 
Dissecting physical models into metric-free components and metric-dependent parts was originally proposed by Tonti, 
[44]. Application of Tonti's ideas for electromagnetism is fully treated by Mattiussi [31 j. A very good introduction 
in the geometric structure of electromagnetism is given by Bossavit, [6, 5|. But based on the analogies described by 
Tonti, the construction advocated by Bossavit have a much wider range of applicability than just electromagnetism. 
Hodge theory of harmonic forms was described by Dodziuk, [ 1 6 1 . Hyman and Scovel, [ 25 1 , derived mimetic operators 
in a finite difference setting which was later generalized by Bochev and Hyman, [ 3 1 . Mimetic finite difference methods 
are described in Brezzi et al. and Steinberg et al., |[T0ll23ll24ll38ll40l . Arnold, Falk and Winther have described an 
extensive framework in a finite element context, |fl~l|2). The geometric ideas underlying this paper are also extensively 
studied by Desbrun et al, [ 14] and Hirani, E2l and DiCarlo et al., ifTSI . For fluid flow calculations, mimetic methods 
were used by Perot, [36, 35] and the importance of preserving physical invariants was illustrated in Perot's review 
paper l34l . Application of these ideas for spectral elements was described in [26, 39] and application of these ideas 
to Stokes flow can be found in [27, 29]. Extension to compatible isogeometric methods can be found in [11] and the 
connection of isogeometric methods with geometry can be found in [20]. The relation between the Hodge matrix and 
mass matrices in finite element methods is discussed by Hiptmair and Tarhasaari et al., |[2Tll43l . 

The outline of this paper is as follows: In Section [2] we introduce the necessary background on vector fields and 
differential forms. In Section[3]the distinction between forms and pseudo-forms is discussed. In Section|4]a discrete 
representation of integrals is given and the analogy with the continuous forms in SectionPJlis presented. In SectionB] 
we describe how we can convert continuous forms in discrete forms and vice versa. In Section|6]we demonstrate how 
the single grid and the dual grid approach can be used to solve the Poisson equation for volume forms. This shows 
how the method works in practice and also shows that both solutions are equivalent. In Section fj\ conclusions are 
drawn and further applications are discussed. 

2. FORMS AND VECTORS 

This section provides an introduction of the basic elements of differential geometry. For more detailed definitions 
the reader is referred to fT8l[T7l . We emphasize here the distinction between vectors and covectors. Just like the 
vector field is an extension of the vector concept to manifolds, the 1-form is an extension of the covector to manifolds. 

As Burke (UJl puts it: 

Were this a mere change in notation, it would make no sense to change things. It is not a mere change 
in notation, however, but a basic change in the fundamental concepts. The new concepts are better for 
unarguable reasons: they [differential forms] correctly represent a larger symmetry group, and therefore 
correctly represent more features of the real world. 

2.1. Tangent vectors and vector fields 

Before introducing differential forms, we need to define vectors or - more precisely - tangent vectors in a domain 
M. In general, Al is a differentiable manifold. Let y(f) be a curve in M parametrized by t e (-e, e), e > 0, with 



p — y(0) e At. The derivative y(0) is a tangent vector at the point p e At. Note that a vector will in general not lie in 
the manifold At itself. So the conventional image of a vector as an arrow connecting two points in the domain At is 
inadequate. For an n-dimensional manifold we can define n curves through the point p e At which produce n linear 
independent vectors. A collection of n linearly independent vectors, e\\ p , . . . , e n \ p , at the point p spans a linear vector 
space denoted by T^At. Any other vector at p can then be written as a linear combination of these basis vectors, i.e. 



Z a\p) e t [ 



where the expansion coefficients a'(p) are associated to the point p and the particular basis e,| p . If we introduce a 
local parametrization of the point p, we can use the coordinate functions x' for the curves which generate a basis at 
p. Such a basis is called a coordinate basis for the tangent space T p Ai and in this case the basis vectors are generally 
denoted by d/dx'\ , or briefly di\ p . For a coordinate basis a vector is represented as 



= J>0>) 



i=i 



dx< 



If we introduce another coordinate system in which to represent p locally, say (y 1 ,... ,y n ), we have 



!>> = 



dy> 



Since any y' - y\x 



, x"), we have 



2A>s 



n n 



'=1 7=1 



3/ d 

dx' dyj 



So b> - (dyj/dx')b'. The fact that a change of coordinates is effectively the application of the chain rule for differen- 
tiation motivates the notation d, for the coordinate basis. 

The construction of a vector at a point p e At can be done for all points in At, such that we smoothly associate 
with each point a vector. This construction generates vector fields. The collection of all tangent spaces is called the 
tangent bundle 

TM := U T p At . 

peM 

2.2. Covectors and l-forms 

With any linear vector space V we can associate the dual space V* of linear functionals acting on V, i.e. 



Va e V , 



V ■ 



So with T p Ai, we can associate the dual space T* At of linear functionals acting on vector v e T p M. 

a : T p M i-» R . 
a (av + bu) - aa(v) + ba(u) , Vw, v e T^At and a e T* At . 



The dual space itself becomes a linear vector space if we set 

(aa + bB) (v) = aa (v) + bB (v) , Vw, v e T p M and a,/3 e r* M . 

The elements a-,/? e T* M are called covectors at the point p e M. The dimension of the cotangent space dm\T* M - 
dimTpM - n. 

Let e\\ p , . . . , e„L be a basis for the tangent space T p M, then a canonical basis for the cotangent space T*M is 
given by e 1 , . . . , e"| p , where the basis covectors satisfy 

[0 if/*; 

In case the vectors are represented in a coordinate basis, <9,, the basis covectors are denoted by dx J . In this case 
any covector at the point p can be represented as 

a - a\(p) djc 1 ] + . . . + a„(p) dx"\ p . 



Application of such a covector to a tangent vector at p then yields 

n ( n 



o 



dxi 



V7=l 



i= 1 j= 1 ^ 

n 

= Y.aiipyip) . (2.1) 

i=i 

Remark 1. In a given basis, application of a covector to a vector resembles the inner product in Euclidean space. 
However, an inner product depends on the metric tensor, whereas a(v) is independent of the metric. Furthermore, 
an inner product is a bilinear form on a single vector space V, whereas application of a covector to a vector is an 
operation between two distinct spaces: V* X V — > R. That these spaces are really different can be seen when we apply 
a change of coordinates. In this case vectors and covectors transform differently and should therefore be treated as 
different mathematical entities. 

Remark 2. A pictorial representation of a vector is usually in terms of an arrow. Covectors can be represented 
by (n — l)-dimensional hyper-surfaces, (ordinary surfaces in 3D). Duality pairing between a covector and a vector 
then yields the number of times the arrow (vector) pierces the surfaces (covectors). For examples of these graphical 
representations, see Bossavit, [6], Burke, M2V and Misner, Thome and Wheeler, H32V . The ultimate aim of this 
distinctive representation is to emphasize the difference between forms and vectors. 

Remark 3. The cotangent space T* M is isomorphic to the tangent space T p M. But there is no canonical isomor- 
phism which associates elements a e T*M to elements v G T p Ai. It is physics which provides a unique connection 
between vectors and covectors. 

Remark 4. We cannot 'see' a form, i.e. a linear functional, directly, but we can only assess its action on the elements 
of the primal space. The only knowledge we can obtain of a form a is by applying a to various elements v of the vector 
space. A similar thing happens for some physical variables: Nobody has seen force', but we only see the action of a 
force. We use scales and spring balances to measure force. The deflection of the hands on the scale or the extension 
of the spring are used to measure force. So we only have access to the work performed by the force, F(v). It therefore 
seems obvious to represent certain physical variables by forms. 



Just as we extended the tangent space in a point p e At, to all points in At to form the tangent bundle, we can also 
consider the collection of cotangent spaces for all points in At, This defines the cotangent bundle 

r At := (J t;m . 

peM 

An element (section) of T*M is then written as 

n 

a ( ^ = Ya k {x l ,...,x n )dx k . 

k=\ 

Such an element from the cotangent bundle is called a differentiable l-form or a 1-form. The space of differentiable 
1-forms is also denoted by A (At). 

Application of a 1-form a^ to a vector field v assigns to every point in A( a real number, i.e. a (I) (v) is a real- 
valued function on At. We will denote the space of real-valued functions on At by A°(At), the space of 0-forms on 
At. 

Example 1. Consider a 1-form cp^ix 1 ,x 2 ) — tp\(x l , x 2 )dx l + <f2(x l , x 2 )dx 2 and a curve y(s) = ly 1 (s),y 2 (s)), such 
that: 

y : [0, 1] i-> At 

The tangent vectors along the curve are given, as usual, by: 

dy 1 Ay 2 

g{s) = -5 — ^i + — — 5 2 
as as 

The action of iff l ^ on g is then given by: 

n A k 

<P m (.g) = y.<Pk-r- = w(s) 

ti ds 

Where w(s) is a function of s and whose Lebesgue integral is: 



f 

Jo 



w(s) = W 



Remark 5. Note that the less usual notation for Lebesgue integral was used. Typically the notation for Lebesgue 
integral shows the measure term, in this case ds, but that is optional. The Lebesgue measure was omitted since, by 
historical misfortune, ds, the Lebesgue measure has the same notation as a 1-form basis dx l , nevertheless both are 
distinct mathematical objects. This overload of notation is the root of a common confusion between differential and 
differential forms (afunctional) and measures: A form changes sign under a change of orientation, whereas a measure 
remains positive. 

Many integrals in physics are still considered as Riemann integrals, where the ds is interpreted as the limit As — > 
ds. But the Riemann integral is defined only for a limited class of functions and this set of functions is not closed under 
the operation of taking point-wise limits of sequences of functions in this class, H30V . All integrals will be Lebesgue 
integrals in this paper. 

Remark 6. As for the physical meaning of this example, consider y> ( ' as a force. Hence, W is the work done by this 
force along the path y(s) and w(s) is a density of work along the path. The integral along a curve is a metric-free 
operation when the force is represented as a differential form. If the force were represented as a vector, as is the case 
in many elementary physics books, work would be a metric-dependent concept. 

Remark 7. Duality between vectors and covectors is reflexive and therefore we can interpret ExampleU\as a recipe 
to integrate the l-form tp"> along the path y, or to integrate the path y against the l-form iff-', H42V . 



Now that we have defined the 0-forms on A1 and the 1 -forms on At, we can proceed in different ways to define 
£-forms for k > 1 . We can either define exterior powers of a vector space to produce £- vectors and then define £-forms 
as the elements of the dual space of the space of &- vectors. This approach is described in H4l [171 [131 . Alternatively, 
we can define a fc-form as an alternating A:-tensor on the tangent space which maps into the real numbers, HTl! 



a {k) : T P M x • ■ • x T P M 



k copies 

with 

a {k \. . . , v t , . . . , v h ...) = -a (k) (. . . , vj, . . . , Vi, ...), v,- e T P M , i=l,...,k. 

Here we can define the exterior product or wedge product to inductively construct £-forms, for 1 < k < n. Let 
A k {M) and A'(A() the space of £-forms and Z-forms, respectively, with k + I < n, then the wedge product, A, is a 
mapping: 

A : A* (At) x A' (At) -> A A+ ' (At) , k + l<n 

that satisfies the following properties: 

(a w + p (l) ) A y (m) = a (i) A r (m) + p (l> A y (m) (Distributivity) (2.2a) 
(a (k) A/3 (l> ) A r (m) = a ik) A 06 (/) A y (m) ) 

_ a W A yjO A y (m) (Associativity) (2.2b) 

fl ff (i) A ^ (0 = ff (/:) A fl^ (0 = fl(Q' ( * , A j6 <0 ) (Multiplication by functions) (2.2c) 

a (k) A p(D = (_if0J) A ff (« (Skew symmeti-y) (2.2d) 

where a (t) 6 A*(At),/3 (/) e A'(At) and y (m) e A'" (At). 

Consider a sufficiently smooth bounded 71-dimensional oriented manifold M c R" with boundary dM and a 
coordinate system where points in the manifold, x € M, are represented by an n tuple x := (x 1 ,- ■■ ,x"). Let a* 
denote a differential £-form, k < n and A* (Ai) denote the space of differential k-forms or k-forms, then a (i) e A* (A1), 
can be written as 

a w = ^ a i (x) t/x'' 1 A dx h A ■ • ■ A rfx", (2.3) 

/ 

where I - i\,- ■ ■ , /^, and \ <i\ <■■■ <ik<n and where a/ (x) are continuously differentiable scalar functions. 
2.3. Integration of differential fo rms 



A Physical system generally preserves integral quantities and these quantities are topological in the sense that 
they do not depend on a particular coordinate system or the metric employed. The integration of differential 
forms is independent of any metric notions in space, |[T8l Ch.3]; no Riemann metric, dot products or k- 
dimensional volume elements are required. If we would have used vectors for integration all these metric 
concepts need to be taken into account. Integrals are of paramount importance in mimetic discretizations, 
because we cannot represent differential forms and the operations on and between differential forms in a 
finite dimensional setting. But we can represent integrals of differential forms and the integral relations 
between them exactly. 



In R 3 there are four types of differential forms, that is 0-,l-, 2- and 3-forms, as many as the different kinds of 
submanifolds: i.e. points, lines, surfaces and volumes. Differential fc-forms naturally integrate over fc-dimensional 
manifolds. Let a {k) e A k (At) and At* c R", with k - dim(Mk), 



{a {k \M k ):= f a w , (2.4) 

•JM t 



this represents a metric-free duality pairing. 

Differential forms live on manifolds and can be easily transformed from one manifold to another, under the action 
of a special mapping called the pullback. Let <£ : At re / — > At be a mapping between two manifolds, the pullback 
operator, O* : A k (M) — > A* 1 (At re /J is such that, 

f a®= f <D*a». (2.5) 

Jo(A1„/) Jm«/ 

The actual evaluation of the integrals is analogous to the case discussed in Example [T] Choose a parameterization 
of the manifold Mk 

r k : [0, 1]* -> M k , 



and then 

f cP>= f o®= f r*(a (t) ), 

Jm- ^ti([o,i]*) J[o,i]» 

where the integral on the right denotes the usual integral in R*. For more complex manifolds one might use a collection 
of parameterizations of the manifold, ri, for the representation of the manifold in which case the integral is evaluated 
as 

f a»> = V f a W = Vf (r<) *(««) . (2.6) 

Jm* ^ Jt^([0,1F) V J[0,1]< 

The inner-product of £-forms, (■, ■), is a symmetric, bi-linear mapping: 

(•, •) : A* (At) x A* (At) -> A°(At) . 
Let a (i) and/?'*' be two A; forms written in the form ( |2.3[ ), then the point- wise inner product is given by, IfTSl fT71 [131 

{a^,^>)^g i -K..g i ^a iu ... Jt /3j u ...j t , (2.7) 

where g k ' 1 : = (dx k , dx l ) are the components of the metric tensor. 

A combination between the wedge operator and the inner product induces the Hodge-* operator, *, which is a 
mapping: 

* : A* (At) -> A""* (At) (2.8) 

such that: 

a*>A*0®:=(a<»£<*>)<7W (2.9) 

where <x (n) € A" (At) is the volume form, that is the n-differential form defined on an n-manifold M such that J cr^ 
is the volume of the manifold M and we have • 1 = <x (n) . The important role of the Hodge-* with regard to orientation 
will be explained in the next section. 

The L 2 inner product, (■, Ol^aq* is defined as a mapping: 

( v )L*(Ai):A*(At)xA*(At)^R, 

such that: 

v w A*/? (/) . (2.10) 



(^•nw-L^K-X/ 



A crucial operator in differential geometry is the exterior derivative, d, defined as a mapping 

d : A* (At) -> K k+l (At) 



The next equation is the most important equation in mimetic methods and can therefore be referred to as the 
Mother of all equations. It relates the geometric boundary operator to the exterior derivative. 

The generalized Stokes Theorem 

f da (k) = f a (k) «* (da < - k) ,M M ) = (a < - k \8M M ), (2.11) 

where Aik+\ is a (k + l)-dimensional manifold and dAik+i is its boundary, a ^-dimensional manifold, and a* 6 
A k (AT). This represents the generalized Stokes theorem that encapsulates the Newton-Leibniz, Stokes and Gauss 
theorem and hence generalizes the vector calculus V, Vx and V- operators to arbitrary manifolds. Moreover, this 
operator is independent of any metric or coordinate system and satisfies the Leibniz rule d (a® A /3® J = da (i) A y3® + 
(-l)*a (i) A 48®, and is nilpotent, dda* := 0. The pullback operator and the exterior derivative have the following 
commuting property, 

A 4 (W() — ^— ► A k+1 (M) 



1** 
A k (M ref )^A k+1 (M ref ) 



d 



Finally, the last piece of machinery is the co-differential, d*, which is a mapping: 

d* : A* (M) -> A*" 1 (At) , 
that satisfies the identity: 

(aP-V.d'fP) =(dcfl* 1 \fP>) - f tra^Atr*^, (2.12) 

J 3 At 

where tr is the trace operator, which is the pullback of the immersion of boundary of the manifold into the manifold, 
C, with ( : 8M i-> A4. When [ tra (t_1) A tr -kff k) - then d* is the formal Hilbert adjoint of the exterior derivative 
induced by the L 2 inner product ( |2.10[ ): 

(dof- k - l \^) M = (a*" 1 ), d^)^ , (2. 13) 

Opposite to the exterior derivative, the codifferential is a metric-dependent operator, given by, d* = (_i)"(* +1 ) +1 * j* 
The Laplace operator, A, is a mapping: 

A : A* (M) -> A* (M) 

which, in terms of the exterior derivative and the co-differential, is given by 

A = d*d + dd*. (2.14) 

For a 0-form a (0) e A (At) the Laplace operator is given by Aa (0) = d*da (0) . For ra-forms <jj (n) e A" (M) the 
Laplace operator is Aw (n) = dd*<j/ M> . 

3. ORIENTATION 

Now that we have introduced sub-manifolds and differential forms as metric-free ingredients in the evaluation 
of integrals, we need to discuss orientation of space and its geometric objects. Orientation is something we need to 
introduce to do our calculations, but physics should - in a sense - be independent of the choice of our orientation. 

S 



From our earliest physics courses orientation has played a role: the right-hand rule in electromagnetism, the 
outward unit normal, calculation of circulation counter-clockwise, vorticity in aerodynamics is oriented clockwise, 
definition of the vector product of two vectors, etc. Many sign conventions in physics are also due to a preferred 
orientation. 

The basic question with orientation is: 'What happens in the description of physics if we change the orientation?' 
Of course, it will only affect the description of physics and not physics itself. A few examples might illustrate the 
point. 

Example 2. In ExampleUjwe calculated the line integral along the path y(s) — (~y l (s),y 2 (s)), < s < 1. If the l-form 
integrated along this curve denotes force, then the line integral denotes the amount of work, W, done by the force 
along the curve. 

Now suppose that we change the direction (orientation) in which we traverse the curve, i.e. we now take the curve 
yCO — (y (1 — sXy 2 (l _ S))> < s < 1. Then the value of this integral is —W for a conservative force. This is as 
expected. So changing the orientation in the calculation of work for a conservative force leaves the differential form 
unchanged and reverses the integral: Wab = — Wba> if the point A and B are the endpoints of the curve on the manifold 
M. 

Is this generally true? If we change the orientation of the domain of integration, does the integral quantity always 
change sign? The answer is no. A simple counter example is the 3-form mass density, p (3) . 

If we integrate p^ over a positively oriented volume (right-hand rule), we obtain the mass in that volume 



-f 

JMi 



M = p {i) > . 

If we change the orientation of the volume and the differential form p (3) would remain unchanged, the integral would 
be 



M 



f pa = - r p<» < o 

Jm-, JMt 



So a change in orientation would lead to negative mass, which is physically unacceptable. Not because we think of 
mass as being positive - hypothetically mass can be negative -, but because its sign depends on the orientation. In 
order to restore our faith in physics and disappoint those who suffer from overweight, the mass density 3-form should 
change sign when the orientation is changed. In this case, mass will always be positive, regardless of the sense of 
orientation. 

So we have (at least) two classes of differential forms: those that do not change sign when the orientation is 
reversed and those that do change sign when the sense orientation is changed. 

3.1. Orienting geometric objects 

Orientation is an indispensable part the description of physics and physics compatible numerical techniques should 
take orientation into account; either explicitly by assigning an orientation to the computational mesh, or implicitly by 
following sign conventions. Orientation is discussed in any textbook on differential geometry such as lfj"8l §2.8]. A 
very good introduction is also given by Bossavit, Q. 

As we will show in the remainder of the paper, orientation is also an important aspect of computational physics. 

3.1.1. Orientation of a vector space 

Let V be a vector space and let e = (ei , . . . , e„) and f = (f j , . . . , f„) be two bases. Then there exists a non-singular 
transition matrix P such that e = Pf. If det(P) is positive, then both bases have the same orientation and if det(P) is 
negative both bases have opposite orientation. This divides all bases into two equivalence classes. All bases in a class 
have the same orientation. This equivalence relation is reflexive, transitive and symmetric. 

We orient a vector space by selecting a basis from one of the two equivalence classes and call this the positive 
orientation or the default orientation. There is nothing intrinsically 'positive' about this basis, it is just a choice. If 
we take any basis from the same equivalence class, this basis will have the same orientation. If we take a basis from 
the other class it will have the opposite orientation. Denote the equivalence class from which we pick out default 
orientation Or and the other class -Or. Bossavit calls the orientation from the default orientation direct frames and 
the bases from -Or skew frames. 



3.1.2. Orientation of a geometric object 

Consider the basis geometric objects in 3D, the point, curve, surface and volume. We can consider these objects 
as sub manifolds in 3D and therefore they possess charts from the objects to R*, k — 0, 1,2,3. Then we set up 
a coordinate basis in each chart. All coordinate systems within a chart have the same orientation. At the overlap 
between charts, we impose that the Jacobian is positive thus ensuring that in a neighboring patch the orientation is in 
the same equivalence class. If we can do this consistently for the whole manifold, we call the manifold orientable. 
Not all manifolds are orientable. A notable example is the Mobius strip and the projective plane RP 2 . 

3.2. Forms and pseudo-forms 

In the examples above, work and mass, we saw that some variables change sign when the orientation is reversed, 
while others do not change sign. The forms which do not change sign are called genuine forms, while forms which 
do change sign are called pseudo-forms. This is an important distinction between physical variables. We can never 
equate a form to a pseudo form, because a change of orientation would falsify the relation. Before we can equate 
a form to a pseudo form, either the form needs to be 'translated' to a pseudo-form, or the pseudo-form needs to be 
converted to a form. 

The notions form and pseudo-form are taken from differential geometry and can be found in almost any book on 
the subject. In this paper we loosely follow Frankel, [f8] Ch.2] and Bossavit, [6|. One may find other classifications 
in literature which mean exactly the same thing: Tonti, [44|, describes the two classes as configuration variables and 
source variables. In earlier work, 11261 we referred to inner- oriented variables and outer-oriented variables. The 
genuine forms were referred to as inner-oriented variables and the pseudo-forms were called outer-oriented. Bossavit, 
IS) talks about axial vectors or twisted vectors and true vectors. Whatever the names one uses to refer to the two 
distinct types of forms, it is important to acknowledge this difference. Not only for the sake of a sound physical 
description, but also for a proper numerical model. 



The operator which converts forms into pseudo-forms and vice versa is the Hodge-* operator, 1T81 §14.1], 
defined in d2~8ll and dZ9li. 



4. DISCRETE REPRESENTATION OF INTGRALS 

Consider a three dimensional domain M and its associated grid, as shown in Figure [T] The grid not only con- 

Cell complex 








Figure 1: Three dimensional domain (left) partitioned into a collection of points, lines segments, surfaces and volumes (right). 

sists of points as is common in many numerical methods, but also the line segments connecting the points, surfaces 
bounded by these line segments and volumes bounded by these surfaces. The partitioning of the domain (manifold) 
in these geometric building blocks is an instantiation of a cell complex. Loosely speaking a cell complex consists of 
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a collection of sets, Q, in this case the ^-dimensional objects in the grid, and an operator d which maps elements of 
the set Ck into the set Q_i, such that ddCk = 0. For a grid to be considered as a cell complex we impose that if a 
^-dimensional geometric object is in the grid, its boundary is also in the grid, where the operator d is the boundary 
operator. 

If we endow all geometric objects (sub-manifolds) with a default orientation, then we call the cell complex an 
oriented cell complex. 

Let M be a manifold covered by an oriented cell complex, D, or an oriented grid consisting of points, lines, 
surfaces and volumes. We will call the individual A:-dimensional geometric objects in D k-cells r^y, k - 0, • ■ • , n, 
where k denotes the dimension of the object and i is a label to distinguish the different points and lines, etc. Given the 
cell complex D, the space of fc-chains of D, C* (D), is the collection of weighted £-cells. A £-chain, C(^ e Ck (D) is a 
formal sum of £-cells, T(k),i e D, 



cm 



Yj c ' T (k)J- (4-1) 



By formal sum we mean a collection of cells and weight {T^y, c'}, where the c' denote the weights. The addition of 
two geometric objects, say points, is not defined, therefore this weighted collection, the formal sum, should not be 
confused with ordinary summation. 

The boundary operator on A:-chains, d : Ck (D) — > Ck-i (D), is an homomorphism, 

dc m = 8 ^ cW m := ^ Jd (r W;/ ) . (4.2) 

i i 

The boundary of a £-cell, T{k\i is a (A; - l)-chain formed by the oriented faces of r^y. The coefficients of this 
(k - l)-chain associated to each of the faces is given by the orientations. 

dr m = ^ efa-Dj, (4.3) 

,/' 

with 

e\-\ , if the orientation of t^-\)j equals the default orientation 
e! = -1 , if the orientation of t^-\)j is opposite the default orientation 
ej — , if T(k-\),j is not a face of T(k),j 

Once basis £-cells have been chosen, the space of fc-chains, Ck(D), can be represented by a column vector con- 
taining only the coefficients, c', of the chain. That is, there is an isomorphism if/: 

\jj : Ck(D) ^>R P , p = rank(Q(Z))) , (4.4) 

defined by 



ifr(c (k) ) = ifr 






[c 1 ...c p ] r , p = rank(C*(Z>)) , (4.5) 



where the rank of Q(D) is the number of fc-cells in the cell complex D and the c' are the coefficients of the £-chain 
C(£). The £-chain, c^), is printed in boldface, whereas the vector C(#) of coefficients is printed in regular face. 

Remark 8. Instead of taking the individual k-cells, T(k\i> as basis, one may choose to take any other collection of 
linearly independent k-chains, C(k)j, as basis. Both bases are related by T^y = A J C(k),> we have 

c (k) = X C ' T «-'' = Zj c ' X A i 5{k> 'i = Zu ^ g »W . 
' i i 

where c 1 — A J c'. This change of basis allows one to incorporate isogeometric methods in this mimetic framework, 
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Example 3. Consider as an example the 2-dimensional grid shown in Figure p] The points Pi are assumed to be 
sinks. So when something 'enters ' a point we denote this by the positive orientation and when something 'leaves ' the 
point it will have the opposite orientation. The default orientation for the lines and surface is indicated in the figure. 
For the figure we see that the boundary of the line L\ is given by dh\ — Pi - P\. Similarly, we can write down the 




Figure 2: Small oriented 2-dimensional cell complex 
boundaries for all line segments and collect everything in a matrix 



( dLi \ 
dL 2 
dL 3 
dL 4 i 



( -1 1 ^ 
0-11 
-10 10 
-1 1 J 



I Pi \ 
Pi 
Pi 

v Pa j 



<3r (iy = ^ e/r ( o)j 



(4.6) 



Note that we stretch, twist or bend this little grid without changing the connectivity and orientation, this matrix 
relation remains valid. This is a metric-free relation because it is independent of shape and size. 
Similarly, we can relate the boundary of the surface to the line segments 



dS = (1 - 1 -11) 



( U \ 
Li 

u 



f)T, 



(2),/ 



2 e i T ^ 



(4.7) 



Again, deformations of the grid which do not change the connectivity between surfaces and lines or orientation give 
the same matrix representation of the boundary operator. 

For a chain complex we need to have that d o d = and this also follows from the matrix representation of the 
boundary operator given in this example 



ddS = (1 - 1 -11) 



( dLl ) 






r -i 


1 


^ 


( Pl 1 


dL 2 
dL 3 


= d - 


1-11) 




-l 






-1 1 

1 


Pi 
Pi 


{ dL A ) 






l 


-1 


1 j 


\ Pa ) 



= 0. 



(4.8) 



Example [3] shows that on an oriented grid, the boundary operator can be represented by an incidence matrix, 
E(JUt-i), which relates a ^-dimensional geometric object to its (k- l)-dimensional boundary. We will denote the matrix 
representing the boundary operator acting on fe-cells, E^_i). The relation d o d — then reads ^(k,k-\)^(k-\,k-i) - 0. 

Remark 9. The incidence matrices, E^k-i), are a matrix representation of the boundary operator. 
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Remark 10. We have that tlie boundary of the boundary always yields tlie empty set. If the boundary of a geometric 
object is the empty set, this does not necessarily imply that the object itself is a boundary. Only on contractible 
domains we have that when the boundary of a k-dimensional object is empty, that object must be the boundary of a 
(k + 1)- dimensional object. This is the Poincare lemma. 

Dual to the space of A:-chains, Q (D), is the space of £-cochains, C (D), defined as the set of all homomorphisms, 

c (k) : C k (D) -^ R 

(c w ,c w ):=c w (c w ). (4.9) 

Remark 11. Essentially, a k-cochain assigns a physical value to a geometric object. For instance, it can assign mass 
to a 3-dimensional object, flux to a 2-dimensional object, temperature to a point, or work to a curve. In this respect 
duality pairing between a cochain and a chain is the discrete analogue of integration in the continuous setting, \2.4\ 



Remark 12. The value assigned to a k-chain (or k-cell) is assigned to the geometric object as a whole. There is no 
particular point in this object where this value is anchored. So if we assign 1 gram to a 3-dimensional object, 1 gram 
will be the attribute of that object and we do not assign it to the center of mass. The center of mass is a priori unknown 
anyway, because we only know the mass of the object and not its mass density distribution. 



With the duality pairing between cochains and chains, one can define the formal adjoint of the boundary 
operator, the coboundary operator, 6 : C k (D) — » C k+l (D), 

(6c (k \c ik+l) ):=(c (k \dc (k+l) ). (4.10) 

This equation is the discrete analogue of the generalized Stokes equation, (|2.11|. 



Remark 13. The coboundary operator, 6, is the discrete analogue of the exterior derivative, d. Derivatives require a 
certain smoothness on the objects to be differentiated. Such smoothness assumptions are not available in a discrete, 
topological setting. By defining the discrete derivative in terms of the boundary operator - which is well-defined for 
discrete objects - we circumvent any smoothness requirements. 

Just like the exterior derivative, the coboundary operator is nilpotent 56c ik ^ - for all c ( *' e C k (D). This follows 
directly from the definition and the fact that d o d = 0. 

Let C k {D) be the space of fc-chains with basis It^jJ, then a dual basis, It (<:) ''J of C k (D) is given, such that 
t^' ! [T(k),j) - <?j- All £-cochains can be represented as linear combinations of these basis elements, 



C W 



= £c,-T ( * w . (4.11) 



Remark 14. Note also here the resemblance between T®' ! (t^j) = 5'. and dx'(d/dx-i) = 5'- for the duality pairing 
between vectors and forms. The k-cochains therefore acts as the discrete differential k-form. 

Once a basis for ^-cochains has been chosen, any ^-cochain is completely determined by the expansion coeffi- 



cients, c,, in (4. 1 1 1. That is, there is an isomorphism ip 

i// : C k (D) ^ W , p = rank(C*(D)) = rank(Q(D)) 
defined by 



v ; 



[ci ...c p ], p = rmk(C k (D)). (4.12) 



Note that a £-cochain is represented by a row vector, instead of a column vector for the expansion coefficients for the 
chains. Duality pairing of a cochain and chain in terms of the expansion coefficients then simply reduces to 



(c«,c w ) = j>< 
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Again note the resemblance with the duality paring at the continuous level given by (2.1 1. 

With the isomorphism which identifies the chains and cochains with their expansion coefficients, we also have a 
natural matrix representation for the coboundary operator 



rank(c t (D» 
(c w ,dc (/t+ i)} = ^ a (E (t+U) c') 

rank(c t (D)) 

J] (ciE ( k + i,k))c' 

= (6c (k \c {k+1) ) . 

So, when the row vector c ik) contains the expansion coefficients, c,, for the cochain c ( *\ then the row vector c^E^+\,k) 
contains the expansion coefficients for dc^K 

Remark 15. This relation is one of the remarkable properties of mimetic methods. Bear in mind that the coboundary 
operator encodes the action of the gradient (k — 0), the curl (k — I) and the divergence operator (k — 2) at the discrete 
level. Once a basis has been chosen, the matrix operation which represents this operation is given by the incidence 
matrix of the oriented grid. These relations are exact, no approximations have been performed. So the topology of the 
oriented grid determines the action of the grad, curl and div. 

In Section|3]we made a distinction between true forms and pseudo-forms. The operator which switches between 
the two is the Hodge-* operator which was defined in ( |2.8[ ) and ( |2.9| l. This cannot be accomplished on a single 
cell complex. Consider Example [5] This is a 2-dimensional problem, therefore n — 2. The Hodge-* is a bijection 
between fc-forms and {n - £)-forms. In the discrete setting of Example[3] it needs to be a bijection between £-cochains 
and (n - A:)-cochains. However, in this example we have four 0-cells (points) and only one (n - k) — (2 - 0)-cells 
(surfaces). Besides, not only the number of dual cells differs from the number of associated primal cells, but also the 
kind of integrals we want to represent differs, see Section [3] So there cannot be a bijection between 0-cochains and 
2-cochains. This is not a particular problem for the grid presented in Example [3] but for grids (cell complexes), in 
general. 

In order to incorporate the duality between of k- and (n - £)-forms at the discrete level we construct a dual grid. 
The construction is as follows: With every £-cell in the cell complex, D, we associate an (n - k)-cell in the dual 
complex. An example of such a construction is shown in Figurep] The collection of all dual cells constitutes the dual 
grid, Di, which contains the interior dual cells of the dual cell complex D, see [26] for further details. In Figure [3] 



D, 



4- 



4- 







Figure 3: A cell complex (left) and its associated dual grid (right) 

the pseudo-forms (outer-oriented forms) are represented on the cell complex D on the left and the true forms (inner- 
oriented forms) are represented on the dual grid. Note that the dual grid itself is not a cell complex in general, because 
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for a cell complex we require that if an ^-dimensional element is part of the grid, then so is its boundary. That is not 
the case for the dual grid. There are some line segments for which not all boundary points are part of the grid. We 
will see in the remainder of this paper that this 'missing boundary' for the dual grid leads to so-called ghost points in 
finite volume methods and to boundary integrals in finite element methods. 

5. SWITCHING BETWEEN CONTINUOUS AND DISCRETE 

In the previous section we described integrals and integral relations at the continuous level in terms of differential 
geometry and at the discrete level in terms of chains and cochains. In this section we explain how to convert the 
continuous forms into discrete cochains and vice versa. 

In this paper we focus on the Laplace operator given by (|2.14i. In Section 13] we presented this operator at the 



continuous level in terms of differential forms. In Section HI a discrete description was given. In this section we are 
going to introduce the operations which will allow us to switch between the continuous formulation and the discrete 
formulation. 

5.7. From continuous to discrete 

The reduction operator, % : A k (vVt) — > C k (D) maps differential forms to cochains by 

(*aP\T m ):= f /» = (/»%). (5.1) 

•Jtau 

Then, for all C(j) e Q CD), the reduction of the £-form, a (k} e A* (At), to the £-cochain, a {k) e C k (D) is, 

a*> (c w ) := (*o«\ c w ) = £ J (XaF>, r m ) = £ <? f KeF> = f a®. (5.2) 

The reduction has an important commuting property with respect to differentiation in terms of exterior derivative 
and coboundary operator, 

<Rd = 6K. (5.3) 

Remark 16. Commuting relations are essential to mimetic methods. Essentially they state that operations at the 
continuous level are mimicked by equivalent relations at the discrete level. In this case, it makes no difference whether 
we take the derivative and then convert to discrete variables or first map to discrete variables and then take the 
discrete derivative. 

5.2. From discrete to continuous 

The other crucial operator is the reconstruction map, I : C k (D) i — > A* (M). This operator needs to have the 
following commuting property, 

dl = 16 . (5.4) 

The map I is an injective, but not surjective. 

Remark 17. The reconstruction operator is essentially an interpolation operator which interpolates the discrete 
cochains to continuous k-forms. The commuting property \5.4\ states that we can reconstruct and take the derivative 
or first take the discrete derivative and then reconstruct. It should give the same differential form. This condition 
poses severe restriction on the admissible reconstruction functions. For triangular elements these reconstructions are 
known as Whitney forms, [6]. 

Furthermore, the reconstruction operator must be the right inverse of %, so %I = Id on C k (D) and it should 
approximate the left inverse of %, so IH — Id + (h p ). 

Remark 18. The first condition, 1U = Id, is a consistency condition. The second condition, 1% = Id + 0(h p ), is an 
approximability condition. Note that the O (hP)-term is in the kernel of the reduction operator, i.e. < RO(h p ) = 0. 
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5.3. Mimetic projection 

The mimetic projection, 777,, is given by 



n h :=Io<R. (5.5) 

A*(A0 ^A£(M;Q) 




where Ajj (At; Q) denotes the range of of the reconstruction, I, of ^-cochains associated with space of ^-chains, 



Due to the commuting properties of reduction and reconstruction we have 

n h d = im = 15'R = dJ7? = d;r A . (5.6) 



Remark 19. 77z/.s commuting property is the most important one of all and it has important consequences. For 
instance, let a® be in the null space ofd, i.e. da^ = 0*- \ then its projected form, a, , will also be in the null space 
of d. Conforming null spaces play an important role in mixed finite element methods, l\29\ l9l/, and the existence of 
discrete potentials, [38 J. 



Another important property is the commuting relation between the mimetic projection, 717, and the pullback 
of a map <t>* 

n h <b* = 0*7r ft . (5.7) 

This will allow us to do all computations on a reference element instead of working directly in curvilinear 
coordinates. For a proof of this property, see 



Having defined the reduction and reconstruction operators for the cell complex a similar derivation can be done 
for the dual cell complex, D, see its construction on page 
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The space Ajj (M; Ck) is the space of discrete fc-forms, with £-cochains associated with the dual fc-cells. 

A* (M; C k ) = n h k k (At) := IKA k (At) , (5.8) 

where Q is the space of ^-chains on the dual complex D, %, is the reduction of differential forms on the dual chains 
and I constitutes the reconstruction of differential forms from cochains defined on the dual complex. 
A discrete wedge product is introduced such that A/, : K k h x A^ — > A^ +/ , given by 

afA h ^:=n(a^A^). (5.9) 

Remark 20. The discrete wedge product satisfies all properties of the continuous wedge product except associativity. 
Associativity can be restored by using a projection nj v h < h, i.e. a mimetic projection on a refined cell complex. The 
precise refinement needed to fully represent the wedge product at the finite dimensional level depends on the number 
of terms in the wedge product and the type of reconstruction forms. 

There are essentially two different ways in which we can incorporate the action of the Hodge-* operator in the 



finite dimensional setting. Either we can use the definition of the inner product for differential forms, (2.9 1 or we can 
use reconstruction of the cochains to obtain differential form, then apply the Hodge-* operator and subsequently we 
reduce the result on the topological dual grid, i.e. 'R-k I or 7? • I. In the former approach we make use of 



,L1M J M h 



K^fU= I «fA*<. (5.10) 



In the definition of this inner product, the Hodge is taken of the second argument and therefore the inner product 
implicitly. Methods based on the use of this inner product will be referred to as the single grid method. 

The method where the Hodge is applied to reconstructed cochains and reduced onto the dual grid will be called 
dual grid method, because for this approach an explicit dual grid needs to be defined. 
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5.4. Basis forms 

Let M be the computational domain decomposed into M non-overlapping, possibly curvilinear quadrilateral or 
hexahedral closed sub-domains, Q m , 



M 



M 



(J Q m , Qm n Qi = <9g m n dQ h m * I, 



(5.11) 



m=\ 



where in each sub-domain a Gauss-Lobatto mesh is constructed. The collection of Gauss-Lobatto grids in all 
elements constitutes the cell complex D. For each element Q„, there exists a sub cell complex, D m . Note that D m n D\, 
m + I, is not an empty set in the case they are neighboring elements, but contains all fc-cells, k < n, of the common 
boundary. 

Each sub-domain is a map from the reference element Q re f - [— 1, 1]", n - dim(AI) using the mapping O m : 
Qref — * Qm- All differential forms defined on M m are pulled back onto the reference element using the following 
pullback operation <D* : A* (Q m , C k ) -* A* (g re/ , C k ). 

The cochains are approximated using piecewise polynomial expansions on the quadrilateral or hexahedral ele- 
ments using tensor products. Thus, it suffices to derive the basis forms in one dimension and afterwards construct 
the n-dimensional basis forms. Furthermore, because of the commutation between the projection operator and the 



pullback, (5.7 1, only the interpolation for the reference element is shown. 



Consider a 0-form, a (u> e A (Q re f), where Q re f := £ e [-1, 1], on which a cell complex D consists on iV + 1 
nodes £,, where -1 < £q < • • • < (jtf < 1, and .A edges, T(iy = [£j_i,£], of which the nodes are the boundaries. 
Corresponding to this set of nodes (0-chain) there exists a projection,^;,, using the .A'' 1 order Lagrange polynomials, 
hi (£), to approximate a 0-form, 



n h a ( - 0) (£) = Y J a i h i (f). 



/=() 



(5.12) 



o 





GLL nodal interpolation 
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Figure 4: The Lagrange polynomials, hi(g), for JV = 4 

Lagrange polynomials interpolate nodal values. Thus, they are suitable to reconstruct the cochain a (0) = Ra^. 
These polynomials are constructed such that their value is one in the corresponding point and zero in all other mesh 
points, see Figure [4] 



^©=/^)= J J *;=" 



i + p 



(5.13) 
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Similarly for the projection of 1 -forms Gerritsma [19| and Robidoux [37] derived 1 -form polynomials called edge 
polynomials, e, (£) e A^ \Qref\ C\\ More details can be found in 



!-l 



e,-(0 = «j(0#, with e, = -^ 



k=0 



dh k 
~dj' 



(5.14) 



The cochain corresponding to the line segment (1-cell), T(i),j is given by m,- = a, - a,-] and so m (1) = r5a (0) is the discrete 
derivative operator in ID. This is a purely topological operation and dla^ — I6a^°\ The 1-form edge polynomial 
can be separated into its polynomial and its basis, 



GLL edge interpolation 






-0.5 



0.5 



Figure 5: Edge polynomial 1-forms, e,-(f), for N = 4 

The edge functions are constructed such that when integrating e, (£) over a line segment it gives one for the 
corresponding element and zero for any other line segment, see Figure [5] 



Ke t (0 






(£> 



1 if i - p 
if i + p 



(5.15) 



The full connection between continuous and discrete counterparts can be summarized in the following diagram, 

C°(D)— 6 ' ~ ° ' ~ ° 

I°y 



A (At) 

A 3 (AT) 



£0 



K 3 



A 1 (At) 

* 

I- 

A Z (M) 



C l {D) 

n 1 ' 



*/, 



C 2 (D) *- 
ft 2 



r 

A 2 (A1) 
• 

6 | P 
A 1 (AT) 



C 2 (D) 



II 



C 3 (D) 



C 1 (D) <- 
7? 1 



A 3 (A1) 

* 

r 

A u (At) 



C°(D) 



The mimetic framework uses Lagrange, /?,(£) e //A°(A( re f), and edge functions, <?,(£) e L 2 A'(A( re f), for the 
reconstruction, I, where the latter is constructed using the former; i.e., from the finite dimensional 0-form ^/,a (0) = 
Zlo aMO e A°(Al re f ; C ), we define 7r,^ (1) e A*(AW; CO, such that 



^ (1) = ^da (0) = £v,^), 



where b, - a, ;— a,-_i . Because we consider tensor products to construct higher-dimensional interpolation, it is sufficient 
to show that the projection operator is bounded in one dimension, ||26ll . 
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(5.16) 



All the interpolation functions were defined in one dimension the extension to multidimensional is straightforward 
by means of tensor products, 

L U,k (6 V, = Ui (£> ® hj (77) ® h k (0 , h (0 ® ej (1/) ® A» (0 , A, © ® A; (/?) ® e* (£)} 
5 St (£ 7. = [^ (0 ® «j fa) ® «* (0 . «» (0 ® ^ (9) ® «* (0 . «< © ® e; (j/) ® A* (£>} 

The approximation spaces are spanned by combinations of Lagrange and edge basis forms given by, 

n I l(\\ \N,N,N 

A° CM») = span {/*?. ) 

1 U n\ \ \N,N,N 11 t-\\ \ \N,N.N 11 ii\ \ \ N.N.N 

A' (M,„) = span \[L m , ) x span (L (I) , ) X span \(L a) , ) 

h m ' F IV Wll<=lJ=0,fc=0 F IV i,J,kJ2)i=0J=\,k=0 F IV i.M/3 J, ■=(),;=(«= 1 ,< i-n 

-i 1/ m\ lAWV 11 nn |AWV 1/ (in \N,N.N W-Ll) 



^o^-^KOiU 



=1,/= u=o 

AL/VQV 



=1,7=U=1 

Using the Lagrange polynomials and the associated edge polynomials on the dual grid and applying a tensor 
product construction, basis forms for the dual complex is constructed similarly. 

6. POISSON EQUATION FOR VOLUME FORMS 

In this section we want to illustrate how both the dual grid and the single grid approach can be used in practice. 
As a test problem we take the Poisson equation for a n-form. The dual grid approach will resemble a staggered finite 
volume method, whereas the single grid approach leads to a stable, well-posed mixed finite element formulation. 



The Laplace operator was defined in (2.14i. Since d(J- n) = for an n-form, the Poisson equation for an n-form 
reads 

Aw'" 1 = f n} <=> ddV" 1 = f in) . 

Here we present results for n = 2 but the method can be readily extended to other values of n. 

6.1. Dual grid approach 

In the dual grid approach we make use of the fact that d* = (-1)"(* +1 ) +1 * d*, which for n = 2 always yields 
d* = - • d*. In the dual grid approach we make use of two dual grids as discussed at the end of Sectionfflon page 14 
We assume that any quadrilateral element, M m , in the (x,y)-plane is obtained from a map <D> : (£, rj) e [-1, l] 2 — > 
(x,y) 6 M m . Then the pullback O* maps forms in physical space, M,„, to forms on the reference element [— 1, l] 2 . It 
therefore suffices to explore the analysis on the reference domain. The 2-forms (J- n > and / (n) will both be represented 
on the primal grid, i.e. on the cell complex shown top left in Figure [6] In terms of the basis functions introduced in 
Section|5] these forms are expanded as 

N N N N 

4 2) = Z Z Oif^M and fit 2 ) = Z Z M-^/W . 

'=1 7=1 '=1 7=1 

where N denotes the number of surfaces in the x- and y-direction in primal grid of Figure [6] Here the coefficients a>y 
are the unknown coefficients and fii is given by 

fa= f r%v 2) . 

Jft-i -''77-1 

If we use these two expression in -d • d • ay, - / (2) , then the first operation we need to apply is the Hodge-* 
operator. This gives 

N N N N N N 

W 2) = • Yj Z ^MOejiTf) = Z Z UiM&ejOi) = Z Z UiMOhiv) ■ 

i=\ j=\ (=1 7=1 i=\ 7=1 
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Figure 6: Primal cell complex (top left) and its associated dual grid (top right). Since the dual grid is not a cell complex, we complement this grid 
with the minimal number of ^-chains to complete the complex. This completion is the boundary of the dual complex (bottom right). The dual 
boundary is itself a complex and its dual is the primal boundary complex (bottom left) 



Here we used ( |5.14| i to write e,(£) = e(£)d£, we use that •d£d77 = 1 and that we can represent this 0-form exactly 
on the dual grid, since we have as many points on the dual grid as surfaces on the primal grid, by construction, see 
Figure [6] Whenever we refer to the dual grid we put twiddles, r , on the coefficients and the basis forms. 

Remark 21. Let the 2- form be given by (J- 2) — o>(£, if) d£d77, then in 2D -kaP' = o(j;,rj), because for (g,rj) e [— 1,1] 2 , 
the metric tensor is the identity. So the function' a>(^, rj) remains exactly tlie same. But only A^&q is chopped away 
from the expression. In that sense, the -k-operator can be seen as an identity operator - because it does not change the 
functions oj(^, rj) - but only associates this function ' with a different geometric object. This is generally not true in 
physical space M m , where the function u)(x,y) will change when the Hodge is applied. The pullback therefore provides 
the metric connection. In the above case, initially the 2-form was associated to surface and after the operation it is 
associated to points. What is less obvious is that the pseudo-form a/ 2 - 1 is converted into a true form tl/ 0) . The same 
happens for the finite dimensional form or, . 



After the *-operation we are on the dual grid and there we want to apply the exterior derivative. This would give 



us 



JV+l N+l 
«"-0 7=0 



N+l N 

'=i ;=i 



N N+l 



&i,j-i)hi(g)ej(i]) 



(6.1) 



^i-ij)ei(^)hj(?]) + ^ ^/fri,j 
i.i j=i 

Note that here we refer to the values a>oj and ci> N+1 j, for j = 1, . . . , A^ and to w,o and 0>yv+i> f° r i — I, ■ ■ ■ ,N. That is, 
we refer to values in points which are not part of the dual grid, see top right plot in Figure|6] It was already remarked, 
that the dual grid is not a complex and we need to add 0-cells and 1 -cells to turn this grid into a cell complex. 
These additional cells are depicted on the bottom right in Figure [6] If these additional points constitute points on 
the boundary of the domain and if Dirichlet boundary conditions are prescribed, then these additional unknowns in 



(6.1 1 have a known value. If not, we simply add these boundary points to the scheme to have a cell complex for the 
dual grid. If the element depicted in Figure|6]is one of the many elements in the spectral element mesh, we use these 
additional points to connect the solution between elements. 
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Remark 22. When homogeneous Dirichlet boundary conditions are employed, the differences (u>ij — c&i-ij) and 
(cjjj — Wij-i) essentially constitute the action of the coboundary on the dual grid as discussed in Section$4\on page 14 
The incidence matrix on the dual grid is connected to the incidence matrix on the primal grid by ^( n -k,n-k-\) — E£ kl) . 
And the incidence matrix on the primal grid was determined by the connectivity and orientation of the primal grid. So 
here we explicitly see that the topological structure of the mesh in fact determines the derivatives on both the primal 
and the dual grid. 

Remark 23. In general, the approximation ofn-forms, a/"', described in the way above, will lead to a discontinuous 
solution between elements. In fact, the trace ofn-forms is not defined. The connectivity is obtained through •k<xr n >. In 
staggered finite volume methods, the introduction of these additional boundary points is also common and they are 
usually referred to as 'ghost points'. In the single grid approach, boundary values also need to be added which leads 
to the boundary integrals in weak formulations, see RemarkW5\ The language differs, but the geometric structure is 
identical. 



Next we apply the •-operator to (6. 1 1 to obtain an expansion of the form 



N N N N 

• d * 4 2) = - 2 Yj iifwwv) + Z Z Ah^w ■ (6 - 2) 



The Hodge expresses the two components in ( 6. 1 1 on the dual grid in exactly the same components on the primal grid. 
So in that sense, the •-operator is just a change of basis. The main thing, however, is that it converts the true form 
(|6.1|l into the pseudo-form (|6.2|i. 



Finally, we take the exterior derivative of ( |6.2| > to obtain 

N N 

d * d * ">? = Z Z («u - «tu + <, - d) e ^ )e ^ ■ 



'=1 7=1 



(2) 



For the Poisson equation for a 2-form, we need to equate this to f h which gives 

N N 

Z Z dj ~ *tu + < ~ <j-i ~ A/) 'tfrjto ■ 

1=1 7=1 

Since the basis functions e,(£) are all linearly independent, the only way in which we can satisfy this equation is by 
setting the expansion coefficients to zero, so 

ij - iUj + ih - «u-i = fu ■ 

Given the explicit location of the dual grids, we can set up a matrix representation of the •-operators in which case 
the Poisson equation for the 2-form on dual grids can be written as 

E (2jl) H u Ef 2I)J H°'V(4 2) ) = Wf ) • ( 6 - 3 > 

Here the incidence matrix Ep.i) is the incidence matrix of the oriented primal cell complex, which depends on the 
connectivity of the grid and is therefore purely geometric. The incidence matrices are independent of the basis 
functions. The Hodge matrices do depend on the basis functions, which is another way of saying that the Hodge is a 



metric-dependent operator. The map <// is the map discussed on page 13 and given by (4.12 1. If we deform the grid, 
the incidence matrices will remain the same and only the metric -dependent part, i.e. the Hodge matrices will change 
as discussed in [7 , 28 1. 

6.2. Single grid approach 

The dual grid approach uses two grids and an explicit construction of the Hodge matrix. In the single grid approach 
the action of the Hodge operator is incorporated implicitly by using an inner product. The exterior derivative can 
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be exactly represented on a grid due to the commutation relation (5.6i. But we do not have such a relation for 
the codifferential opera tor. In the single grid approach, we therefore want to convert codifferentials into exterior 
derivatives using (2.12 1, In order to do so, we need to rewrite the equation dd*w (2) = / (2) in an equivalent first order 
system given by 

dd*w (2) =/ (2) <=> | (6.4) 

ld^=/< 2 > 

If we take the inner product of the first equation with an arbitrary 1-form, v^ 1 ' and the inner product of the second 
equation with an arbitrary 2-form p (2 \ we obtain 



< ' /L 2 A' V ' /L 2 A' 

"' ;i( 2 A = (f (2 \p (2) ) 







(dqvpv>) L 



(6.5) 



>L 2 A< 



Now we can apply (2.12 1 to the second term in the first equation to obtain: Find q m e HA 1 and <J- 2) e L 2 A 2 such that 
Vv (1) e HA 1 and Vp (2) e L 2 A 2 we have 



(6.6) 



' (P, vO) L2A , - {a 2 \ dv(») L2A2 + j m tr v(» A tr * ^ = 

. K^ (2) k, = (/ (2 v 2) k. 

This formulation is equivalent to the mixed formulation in finite element methods, [9|. 

Remark 24. Where in the dual grid approach we needed to impose that the forms were locally integrable, we need 
to put more strict conditions on the admissible solution for the single grid approach. The solution must belong to 
suitably chosen Sobolev spaces. 



Remark 25. Dirichlet boundary conditions imposed on a/- 2 ' appear through the boundary integral (weak imposition 
of boundar y con ditions). Here we see that we need to prescribe tr ■*• a/ 2 ', just as in the case for the dual grid approach, 
see Remark 22 If we prescribe tr -k w <2) , the boundary integral can be transferred to the righthand side. If Neumann 
conditions are prescribed, i.e. trq^ v> is given, then trv { • ' — and the contribution of the boundary integral vanishes. 



Well-posedness of the weak formulation can be found in any book on mixed finite element methods, for instance, 
@. Since the exterior derivative maps A 1 onto A 2 , we can always find a particular q f such that dq f = / (2) . Let 
Z - {v (1) e HA 1 1 dv (1) = }, then if we restrict the solution space and the space of weight forms to Z we obtain: Find 
<7 (1) € Z, such that for all v (1) e Z, we have 



JdM 



The bilinear form (g (1) , v (1 M 2 is trivially coercive and bounded. Consider the space Z 1 ' - {« (1) e ^A 1 1 (a (1 \ v {l) ) L 2 M , Vv (1 ^ e 
Z )■ Since the map d : Z 1 ' — > L 2 A 2 is a bijection for n - 2, we have the Poincare inequality that for every p {1) e L 2 A 2 
there exists a unique Vp e Z 1 ' c HA 1 such that dv p = p (2) and there exists a constant c p > such that, |[T1 l2l l29l 



(2),| 



v^W < c p \\p w \\ L 2 A 



Then it follows that, E9ll4l 



sup 

v<»eL 2 A 2 



f/ 2 Uv (1) ) , , ( P {2 Kdv ( ' } ), 

V ' /L 2 A 2 _ V ' P )L 2 



L 2 A 2 



l/(l)|| 



llv?W 



"^6- > l||^n 2 
c 



IhtfW 



L 2 A 2 



'P 



Therefore, the inf-sup condition, (5), is satisfied at the continuous level. 
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If we restrict the infinite dimensional function spaces to finite dimensional conforming subspaces using the 
mimetic projection n%, we obtain the discrete equation: Find q h € HK x h and to. e L 2 A^ such that Vv^ e //A ; ' 
and Vp h e L K h we have 

(C ^ AI - (*?>. dvl») L2A2 + j m tr v« A tr * 4 2 > = 
Well-posedness of this discrete formulation follows directly from the fact that Hh l h c H A , i 2 A^ c L 2 A 2 and 



2/i = {vjj e //A ; ' | dvj t -0]cZ The latter property is a direct consequence of (|5.6 1 



If in 


the continuous setting, 


the 


equation dg (1) = 


/ (2) possesses a 


solution g (1) , 


then the discrete solution 


< 


= /^ has a solution ni,q m , 


m 



















= 717, 


(/ (2) - 


-6q^ 


^/ (2) - 


- cW 1 ' 


_ f (2) 


"< 





Remark 26. The inf-sup condition is not really a 'condition', because it is identically satisfied by the mimetic projec- 
tion. For Stokes flow this is proven in H29V . 



The discrete matrix equation corresponding to (pTTb has the following form 



M 2 E 



M 1 

(2,1) 










tr • a> 



(2) 



W) 



UK 



(6.8) 



in which M l is the mass matrix associated with the 1-form basis forms, ( 5.16 l and M 2 is the mass matrix for the 
2-forms, (5.16 1. The incidence matrix, E(2,i), denotes differentiation at cochain level and is again determined by the 



geometry. The map if/ is the isomorphism which connects the cochain to its expansion coefficients, (14. 12b. Elimination 
of \jtqi ) from this system yields 

-l 



M 2 E (2>1) (M'Y E (2>1) (M 2 )iA(4 2) ) = M 2 ^) , 



or after pre-multiplication by the inverse of M 2 

E (2il) (m 1 ) -1 E (2 ,d (m 2 )«a(4 2) ) = fUf) 



(6.9) 



If we compare this with the discrete equation for the dual grid method, (6.3 I, we note that in the single grid method the 



mass matrices and their inverses play the role of the Hodge matrix. This relation between mass matrices and Hodge 
matrices was also observed by l5l l3l IT1 I2TI l43l . 

6.3. Results for the Poisson problem for volume forms 

In this section we will present results for a sample problem using both the dual grid and the single grid. As 
already remarked in Example [3] the incidence matrices remain invariant under distortion of the grid. And (6.3 i and 
(6.9 1 reveal that differentiation is determined by these incidence matrices. When we deform the grid, only the Hodge 



matrices in (6.3 i and the mass matrices in (6.9 1 will change, flT] [33]. Therefore, for this test case we consider three 
meshes obtained by a transformation of the unit domain [-1, l] 2 to curvilinear coordinates. The mapping is given by 
(x, y) - <!>(£, ?]), with 



x(tj, rj) — £; + c sin(7r£) sin(7T77) 
y(</j, rj) = 7] + c sin(7r£) sin(7r?7). 

For c = we obtain the orthogonal Gauss-Lobatto grid. The grid obtained for c = 0.0, c — 0.1 and c 
displayed in Figure [7] 

the exact solution for this problem is given by 

Qf <2) = sin(27rjc) sin(27ry) dxdy . 

For these tests Dirichlet boundary conditions corresponding to the exact solution are prescribed. 



(6.10a) 
(6.10b) 

0.2 are 

(6.11) 
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Figure 7: Deformed geometry mesh for deformation coefficient c — (left), c = 0. 1 (middle) and c = 0.2 (right), for 3x3 elements of order N = 6. 

6.3.1. Results for dual grid approach 

In Figure [8] /z-convergence is shown for polynomial degree N = 1,2,3, respectively, for the meshes shown in 
Figure U\ 
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Figure 8 : Plots of A-convergence of the L 2 error for N = 1 , 2, 3 on the meshes shown in FigureMfor the dual grid approach J6.3J and for the single 
grid approach l |6.8) . Left plot mesh c = 0.0, middle plot c = 0. 1 and right plot for c = 0.2. The straight lines indicate the slope of the expected 
convergence rate. 



The /z-convergence curves for the dual grid approach and the single grid approach are plotted in Figure [8] and 



(2) 



,(!) 



are indistinguishable. The expected rate of convergence for both the solution or. and the fluxes q h are expected 
to be Af + 1, where Af is the polynomial degree, which is confirmed by the results shown in Figure [8] In Figure [9] 
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Figure 9: Plots of p-convergence of the Lr error on a 2 X 2 mesh and a 4 X 4 mesh on the meshes shown in FigureMfor the dual grid approach 
|63land for the single grid approach J6.8) . Left plot mesh c = 0.0, middle plot c = 0.1 and right plot for c = 0.2. The black line is the interpolation 
of the exact solution. 



/^-convergence is plotted for both the dual grid approach and the single grid method. Both methods display identical 
convergence rates. The solid line in FigureMis the L 2 -error of the interpolation of the exact solution. Despite the fact 
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that the convergence rates are exactly the same for both methods, the solutions are not identical. 
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Figure 10: Contour plot of the absolute value of the spatial difference between the dual grid solution and the single grid solution on a 2 X 2 grid, 
with p = 3 (left) and p = 5 (right). 



Figure 10 displays the absolute value of the spatial difference between the solution obtained with a dual grid 



,(1); 



approach and the single grid approach. While the fluxes qi in both methods will be slightly different, the conservation 
law dq h = fl is satisfied for both methods as shown in Figure 11 and Figure 12 For this particular test problem 
_f (2) = -87r 2 sin(27nc) sin(27ry) dxdy. 
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Figure 11: L°°-error for the dual grid approach in the equation dq, — K for c = 0.0 (left), c = 0.1 (middle) and c = 0.2 (right) 





10 15 20 25 



10 15 20 25 




15 20 25 



Figure 12: L°°-error for the single grid approach in the equation &q h — K for c = 0.0 (left), c = 0.1 (middle) and c = 0.2 (right) 
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7. SUMMARY AND OUTLOOK 

In this paper we presented two formulations for the numerical solution of the Laplace equation: One in which the 
action of the Hodge operator is explicitly performed using two topological dual grids and one where the action of the 
Hodge is embedded in the definition of the inner product. In both methods the discrete derivative, which is the formal 
adjoint of the geometric boundary operator, is explicit in terms of the incidence matrix of oriented grid. 

In terms of the L 2 -convergence rate for the Poisson problem for volume forms, for h- and /^-refinement, both 
methods are optimal and the convergence plots show indistinguishable curves on orthogonal and curved meshes. 
The actual solution, however, differ. The major difference between the primal-dual grid formulation and the single 
grid mixed formulation lies in the discretization of the codifferential operator. Despite these difference, conservation 
equation dq h = K is satisfied up to machine precision. 

Current work focuses on eigenvalue problems for the Laplace operator applied to A:-forms and the discretization 
of the Lie derivative in order to model convective behaviour. 
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